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Abstract 



Kinetic Particle In Cell (PIC) methods can extend 
greatly their range of applicability if implicit time 
differencing and spatial adaption are used to address 
the wide range of time and length scales typical of 
plasmas. For implicit differencing, we refer the reader 
to our recent summary of the implicit moment PIC 
method implemented in our CELESTE3D code [G. 
Lapenta, Phys. Plasmas, 13, 055904 (2006)]. In- 
stead, the present document deals with the issue of 
PIC spatial adaptation. Adapting a kinetic PIC code 
requires two tasks: adapting the grid description of 
the fields and moments and adapting the particle de- 
scription of the distribution function. Below we ad- 
dress both issues. First, we describe how grid adap- 
tation can be guided by appropriate measures of the 
local accuracy of the solution. Based on such infor- 
mation, grid adaptation can be obtained by moving 
grid points from regions of lesser interest to regions 
of higher interest or by adding and removing points. 
We discuss both strategies. Second, we describe how 
to adapt the local number of particles to reach the 
required statistical variance in the description of the 
particle population. Finally two typical applications 
of adaptive PIC are shown: coUisionless shocks and 
charging of small bodies immersed in a plasma. 



1 Introduction 

Methods to adapt particle-in-ccU (PIC) kinetic 
plasma calculations are very valuable in the study of 
multiple-length scale problems. Typically, multiple- 
length scale problems present small regions of 
stronger gradients embedded in large systems. In 
such conditions, computational efficiency is achieved 
best by focusing the attention in the regions of inter- 
est. 

In PIC methods it is not sufficient to use adaptive 
grids with finer spacing in the regions of interest, it is 
also necessary to rezone the number of particles. By 
particle rezoning we define the operation of increasing 
the number of particles in regions where higher accu- 
racy is required, and of reducing the number of par- 
ticles where lower accuracy can be tolerated. Finer 
grid spacing leads to a better description of the elec- 
tromagnetic fields, but particle rezoning is needed to 
gain a better description of the plasma dynamics and 
a reduction of noise [?] . 

Particle rezoning can also be beneficial to keep the 
load of work uniform on a per cell basis, a feature of 
crucial interest in a correct load balancing in parallel 
implementations . 

In the present work, we review our work in the field, 
without attempting to present a complete coverage of 
the literature. The paper is organizes as follows. 

Section 2 reports general comments regarding the 
task of PIC adaption, discussing in particular the link 
of particle and grid adaption and the link of both with 



time differencing. Section 3 deals witli grid adapta- 
tion, while Section 4 deals with particle adaptation. 
Section 3 is organized around the two main task of 
grid adaption: where to adapt and how to adapt. We 
answer the question of where to adapt by proposing a 
posteriori measures of the local accuracy. We answer 
the question of how to adapt by examining two dif- 
ferent possibilities: grid motion and grid refinement. 

Finally, Section 5 presents to types of simulations 
where PIC adaption is tested: shocks and dust charg- 
ing. 

2 PIC Code Adaptation 

Multiple scale plasma physics problems present two 
challenging features. First, in any given region of 
the system processes develop at widely different time 
scales. Electrons and ions respond with scales made 
extremely different by their different masses, and a 
host of different intabilities can develop each with its 
own time and length scales. Second, different regions 
of the system can have widely separated spatial and 
temporal scales. For examples, regions of localized 
strong gradients can arise locally, as is the case of 
current sheets in space systems. 

The normal textbook approach to PIC is unsuit- 
able to the conditions described above. The stan- 
dard PIC is based on explicit time differencing and 
is subject to strict stability constraints. The time 
step needs to resolve both light- wave propagation and 
Langmuir wave propagation: 

cAt < Ax (1) 

LUpeAt < 2 (2) 

regardless to our relevance to the scale of interest. 
The grid spacing needs to resolve the electron Debye 
length: 

Ax < qXoe (3) 

to avoid the so-called finite grid instability [?]. For 
this reason explicit methods need to resolve the finest 
scales everywhere. 

When implicit methods are considered [?], the sta- 
bility constraints ([l] [2] |3| are removed and the local 



spacing can be chosen according to the required ac- 
curacy rather than the need to avoid instability. A 
practical condition that ensure good energy conserva- 
tion in an implicit PIC method requires that the av- 
erage electron population does not travel more than 
one cell per time step: 

Vth,eAt < Ax (4) 

This condition can be satisfied if both grid spacing 
and time step (but not the one only) are chosen large, 
to step over the small and fast scales. The scales not 
resolved accurately are not eliminated but rather de- 
formed and the energy present in them is damped like 
in the physical Landau damping but at an accelerated 
rate [?]. 

We refer the reader to our recent review of the 
imphcit PIC algorithm used in the CELESTE code 
for the details of how the implicit moment method 
is derived and used [?, ?]. In the present paper we 
describe, instead, how implicit PIC methods can be 
adapted in space. 

To adapt a PIC code to a local scale length, we 
need to address two issues, how to change the local 
grid resolution and how to change the local statistical 
description of the particle distribution functions. The 
two issues are described in the two sections below. 
The assumption is made that the host PIC method 
be implicit, so that large cells do not lead to the fi- 
nite grid instability described above. Nevertheless, 
some forms of adaptation can also be of relevance to 
explicit methods and some of the methods described 
below can also be used in explicit codes (for example 
the application in section 5.2 is explicit). 

3 Grid Adaptation 

Grid adaptation can be achieved by grid refinement 
(i.e. adding more grid points) in some selected ar- 
eas or by grid motion (i.e. moving grid points to 
regions of interest from regions of lesser interest). In 
the first case, the adaptive mesh refinement (AMR) 
method [?] is obtained. In the second case, the mov- 
ing mesh adaptation (MMA) method [?] is obtained. 
A specific class of MMA algorithms widely used are 
ALE methods [?]. In all cases we need guidance. We 
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need to know what interesting mean. Often, interest 
is defined based on the knowledge of the solution. 

In many plasma physics problems the regions of 
interest are readily identified. For example, in space 
weather simulations localized regions of strong cur- 
rents are site of topological changes and require a 
high resolution while regions of smooth flow can be 
described by coarse meshes. 

However, in other instances it is not obvious what 
regions require increased accuracy. In those cases, 
we need error detectors to tell us where the error 
is larger. Here we describe a specific error detec- 
tor previously applied successfully in plasma physics 
problems: the operator recovery error origin (OREO) 
detector [?]. 

For AMR codes, the OREO detector provides ac- 
curate and automatic determination of where the dis- 
cretization error is being generated. This knowledge 
is directly used by the AMR method to refine or to 
coarsen. 

For MMA codes, the knowledge of the error needs 
to be supplemented by a method to move the grid. 
Given the error what new grid should we use? To 
answer this additional problem typical of the MMA 
method we also present a new technique citelapenta 
based on the Brackbill-Saltzman approach [?]. 

3.1 Automatic Guidance on Resolu- 
tion Requirements 

In a previous paper [?] , we have proposed a new error 
origin detector based on the extension of the gradi- 
ent recovery error estimator [?] . We have named the 
approach operator recovery error origin (OREO) de- 
tector since it extends to any operator the method 
used for the gradient operator by the gradient recov- 
ery error estimator. Below, we summarize briefly the 
procedure involved in its definition and implementa- 
tion. 

For the sake of definiteness, we shall assume a gen- 
eral N-dimensional grid (where one of the dimensions 
could be time) where a vector field v„ is node cen- 
tered. For notation, we label the cells with c and the 
nodes with n, using further the notation n(c) to indi- 
cate the nodes neighboring cell c and c(n) to indicate 
the cells neighboring node n. 



We consider a general multi-dimensional non- linear 
partial differential operator: 

0(g) (5) 

Equation ^ summarizes the most general operator 
acting on a function q(x) defined on the multidimen- 
sional space X. 

Equation ^ is discretized on a grid with N nodes 

On{qi,---,qN) (6) 

From the discretized field g„ and from the dis- 
cretized operator X„ applied to g„ defined only on 
the grid nodes, it is possible to reconstruct two func- 
tions defined everywhere in the continuum space x: 

= En9n'S'(x-x„) 
0(x) =EO„5(x-x„) 

where 5(x — x„) is the b-spline basis function of order 
£ for interpolation. 

The local truncation error is defined as the dif- 
ference between the linear interpolation of the dis- 
cretized operator applied to the discretized field 
Xq{x.) and the exact differential operator applied to 
the linear interpolation of the discretized field g(x): 

e = 0(x) - Og(x) (8) 

The average local truncation error on any given cell 
c is defined as a norm of the error e. The L2 norm is 
often used: 

where Cc is the average local truncation error over cell 
c and Vc is the cell volume. 

3.2 Moving Mesh Grid Adaptation 

We have recently proposed a new approach [?] to vari- 
ational grid adaptation [?] based on the minimiza- 
tion of the local truncation error defined above. The 
method can be constructed starting from the follow- 
ing equidistribution theorem proven in Ref. [?] 
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THEOREM: In a optimal grid, defined as a grid that 
minimizes the local truncation error according to the 
minimzation principle 



d'^x, 



(10) 



the product of the local truncation error in any cell i 
by the cell volume Vi (given by the Jacobian J — ^) 
is constant: 

CiVi = const (11) 

The equidistribution theorem is apphed solving the 
following Euler-Lagrange equations: 



9- 



d 







(12) 



This approach creates a grid where \ei\Vi is con- 
stant. Note that the equations above are identi- 
cal to the equations used by the Brackbill-Saltzman 
variable diffusion method [?]. The primary inno- 
vation is that the monitor function is now directly 
linked with the local truncation error instead of be- 
ing left undefined. In the typical implementations 
of the Brackbill-Saltzman method, the monitor func- 
tion is defined heuristically by the user. The use of 
the OREO detector proposed here results in a more 
accurate scheme [?]. 

We have applied the grid rezoning described above 
to our MMA magnetohydrodynamics (MHD) code 
GRAALE [?] based on the ALE discretization [?]. 
Here we limit the discussion to the classic spherical 
ID implosion test proposed by Noh [?]. An unmagne- 
tized gas with 7 = 5/3 initially has p = 1, e — 10~^ 
and uniform velocity u = —1 (except in the center 
where u{r = 0) = 0). The problem represents a se- 
rious challenge for Lagrangian calculations and the 
solution is known to suffer from serious wall heating 
due to the use of artificial viscosity to capture shocks. 
Note that we are not using artificial heat conduc- 
tion [?] (a tool to mitigate the wall heating problem) 
precisely to highlight the trouble of Lagrangian cal- 
culations for the present case . 

The results of an MMA calculation using the adap- 
tive grid is compared with a reference standard La- 
grangian calculation. Figure [l] shows the density at 
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Figure 1: Noh's spherical benchmark: comparison 
of the density at the end {t = 0.6), for a Lagrangian 
(dashed) and MMA (solid line) calculation. The ex- 
act solution is also shown (dotted line). 



the end of the Lagrangian and MMA calculation. The 
use of adaptive grid results in a much improved so- 
lution. The reason for the improvement is explained 
by the sharper resolution of the shock achieved by 
the adaptation. As noted in the original paper by 
Noh [?] , a sharper resolution of the shock also implies 
a reduction of wall heating, as observed in Fig. [T]for 
the MMA case. 

The use of grid adaptation based on the OREO 
detector results in an automatic method to increase 
the accuracy of the MMA method. 
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3.3 



Adaptive 
(AMR) 



Mesh Refinement 



To investigate the performance of the OREO detec- 
tor in 2D, we have apphed it to results obtained 
with CLAWPACK [?]. CLAWPACK is a pubhcly 
available software based on an AMR solution [?] of 
the conservation laws. We have applied the code to 
the solution of the gas dynamics equations for the 
Colella's wedge problem [?]. A planar M = 10 shock 
is incident on an oblique surface; the angle between 
the shock direction and the surface is 7r/6. The actual 
computed results at time t = 0.2 for a 240x120 grid 
are shown in Fig. |2]where all the expected features [?] 
can be recognized. 

The OREO detector is computed based on the re- 
sults obtained from CLAWPACK using Algorithm 3. 
The detector is shown in Fig. [3] for a simulation with 
a grid 120x60. For comparison we also provide an es- 
timate of the actual error, computed by difference be- 
tween the solution on a 120x60 grid and the more ac- 
curate solution on a 240x120 grid. Clearly the OREO 
detector is successful in detecting all origins of errors. 
The shocks are all captured; the slip surface rolling 
up under the shock is evident. All features are de- 
tected. 

For reference. Fig. [3]-c shows also a similar analy- 
sis conducted on another possible candidate for error 
detection often used in the literature. The detector, 
which we name warp indicator for convenience, mea- 
sure the local error as the variance among the differ- 
ent values obtained at a node when extrapolating the 
internal energy from the four directions (backward 
and forward along x and backward and forward along 
y). The analysis in Fig.[3]-c shows that the two right- 
most planar shocks are captured well, while the top 
and bow shocks are barely visible. All the structure 
inside the rolling up region within the outer shocks 
is lost: no slip surface is measured and the internal 
shock is also lost. In practice the warp indicator is 
often supplemented by other ad hoc detectors to pick 
up all shocks, but still the rolling up region and the 
slip surfaces are often left undetected. 

The OREO detector does not miss any feature and 
can be used reliably alone without any other ad hoc 
detector. 
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Figure 2: 2D Gas Dynamics (Eulerian form)- 
Colella's benchmark on a 240xl20grid. Density, ve- 
locity and internal energy at the end of a Eulerian 
calculation {t — 0.2). Results obtained using CLAW- 
PACK. 
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Figure 3: 2D Gas Dynamics - Colella's benchmark on 
a 120x60 grid at t = 0.2. Comparison of the global 
truncation error (a) with the OREO detector (b), and 
the warp indicator (c). 

4 Particle Adaptation 

Particle rezoning is needed to increase the number of 
particles in regions where high accuracy is required, 
and to reduce the number of particles where lower 
accuracy can be tolerated. The primary effect of in- 
creasing the number of particles is to reduce the vari- 
ance of the statistical description of the distribution 
function. In a PIC simulation this increases the ac- 
curacy defined as typical in MonteCarlo methods, i.e. 
as the variance of the simulation. 

Particle rezoning must be in effect throughout the 
calculation to constantly keep the local required ac- 
curacy. In multiple-length scale problems, the region 
of interest can move, and particle rezoning must fol- 
low the motion to keep the focus where it is needed. 
The approach followed here is to use adaptive grids 
to follow the evolution of the system [?, ?] and parti- 
cle rezoning to keep the number of particles per cell 
constant. This approach leads to finer grid spacing in 
the region of interest and, automatically, to a higher 



density of computational particles in that region. 

The problem of particle rezoning can be formulated 
[?] as the replacement of a set of N particles with po- 
sition Xp, velocity Vp, charge q-p, and mass rUp, with a 
different set of N' particles with position x^/ , veloc- 
ity Vp/, charge Qp', and mass mp>. The criterion for 
replacement is the equivalence between the two sets, 
defined as the requirement that the two sets must 
represent the same physical system, with a different 
accuracy. This generic definition of equivalence be- 
tween two sets is given practical bearing by specifying 
two rules for equivalence. 

Two sets of particles are considered equivalent if 

in 

1. the two sets arc indistinguishable on the basis of 
their contributions to the grid moments; 

2. the two sets of particles sample the same velocity 
distribution function. 

The first criterion concerns the moments of the par- 
ticle distribution used to solve the field equations. 
The moments are defined at the grid points Xg as 

^s = E^('^5-^f)9pF(vp) , (13) 
p 

where S is the assignment function [?, ?]. In general, 
when nonuniform grids are used, x is the natural co- 
ordinate, i.e. the system of coordinates where the 
spacing between consecutive points is uniform and 
unitary in all directions [?]. The function F of the 
particle velocity characterizes the moment. In ex- 
plicit electrostatic codes, only the charge density is 
required: 

Pg = Y,^3i^p'l'^P (14) 
p 

derived from (1) using F(vp) = 1 and using a short 
notation for 5'g(xp) = S'(xg — Xp). Electromagnetic 
and implicit codes [?] require higher order moments 
like the current density 

Jg = ^S'g(xp)gpVp (15) 
p 
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and the pressure tensor 

Iig = J2 Sg{xp)qpVpVp . (16) 

p 

The first criterion requires the two sets of parti- 
cles to give the same moments relevant to the field 
equations. Note that if this criterion is satisfied ex- 
actly total energy and momentum are also automat- 
ically conserved. The second criterion is more diffi- 
cult to apply in a quantitative fashion. In previous 
work [?, ?], it has been proposed to use the test or 
the Kolmogorov and Smirnov test to verify that the 
particle distribution is preserved. In practice, this is 
not easily achieved. 

In fluid PIC codes, the first criterion is the only 
one to be applied, and general schemes for particle 
rezoning can be derived [?]. In kinetic PIC codes, 
the computational particles sample the real plasma 
velocity distribution, and the second criterion must 
also be imposed. In the kinetic case the choices are 
more limited. For this reason, a simpler approach is 
followed [?, ?]. To increase the number of particles 
per cell, a given particle is split in two or more new 
particles displaced in space but all sharing the same 
speed. The weights and displacements can be chosen 
to conserve exactly the grid moments, and the veloc- 
ity distribution is not altered because all the particles 
have the same velocity. 

Another approach can be considered. A particle 
can be split in the velocity space. The daughter parti- 
cles have the same position but different velocity. The 
advantage of this method is that the charge density 
is not affected. However, the higher order moments 
(current density and energy) cannot be all preserved. 
Furthermore, the velocity distribution is altered. 

To decrease the number of particles, the splitting 
operation can be inverted to coalesce two particles 
into one. The difficulty is that, in general, it is im- 
possible to find two particles with the same veloc- 
ity. For this reason, particles with different velocity 
have to be coalesced. To minimize the perturbation 
of the velocity distribution, the particles to be coa- 
lesced must be chosen with similar velocity. An al- 
ternative approach is to coalesce three particles into 
two, which allows one to conserve both energy and 
momentum [?]. 



In the following sections, we will provide the two 
most successful general techniques to adapt the num- 
ber of particles in a cell. We refer the reader to a pre- 
vious technical description of the various alternatives 
and their merits [?] 

4.1 Summary of the Algorithms for 
Pctrticle Rezoning 

In the previous sections, we derived all the required 
blocks to build algorithms to change the number of 
particles in any given cell. Here, we provide a precise 
algorithmic description of the methods to increase 
the number of particles per cell and to decrease the 
number of particles per cell. 
Splitting Algorithm : 

Given a cell g with Np particles in a ID, 2D, or 3D 
system, any chosen particle (labeled a) with charge Qo 
( and mass obtained from the charge-to-mass ratio for 
the species), position (in natural coordinates) and 
velocity Vq can be replaced by N' particles, labeled 
p' = {1, 2 . . . N'}. In ID, N' = 2 and the new prop- 
erties are qp' = qo/2, Xp> = ± ^/Np (where the 
cell size is unitary), Vp' = Vq. In 2D, N' = 4 and 
the new properties are Qpi = qo/^; Xi.2 = Xq ± l/Np, 

X3,4 = Xo, J/1,2 = Vo, 2/3,4 = Vo ± 1/-^?/ Vp' = V^. In 

3D, N' = 6 and the new properties are qp' = qo/Q, 

2:1,2 = Xo± 1/Np, X3,...6 = Xo,- 2/1,2,5,6 = Vo, 2/3,4 = 
yo ± l/iVp/ -Zl,...4 = Zo,- 25,6 = Zo± 1/Np. 

Note that the choice of the particle p = o in the set 

of Np particles in the cell g is free. In the result sec- 
tions, we choose the paxticle with the largest energy: 
rUpVp. Algorithm SI preserves exactly the velocity 
distribution function and grid moments. However, 
for quadratic assignment functions the grid moments 
are only approximately preserved (see Section HI). 
Coalescence Algorithm: 

Given a cell g with Np particles in ID, 2D, or 3D 
systems, choose N = 2 particles p = {1,2} close to 
each other in the phase space. Their properties are Qp, 
Xp, and Vp. The two chosen particles can be replaced 
by one particle (labeled A) with qA = qi + q2, = 
(giXi -\- q2X2)/qA, VA = (giVi -\- q2-V2)/qA- 

Algorithm CI preserves the overall charge and mo- 
mentum and the charge density pg but perturbs the 
velocity distribution. Note that one can choose va 
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to preserve the energy, but it is not possible to pre- 
serve energy and momentum together. The crucial 
point of algorithm CI is to choose two particles close 
in velocity and space. A pair search of the two par- 
ticles closest in velocity is usually too expensive. For 
this reason, we perform a diatomic search that sorts 
the particles into two bins and selects the largest bin. 
The binning is repeated in sequence for each spatial 
direction and component of the velocity. The bin- 
ning is continued until the number of particles in the 
largest bin is small enough to use a pair search. 

5 Examples of Adaptive PIC 
Simulations 

To illustrate the possible applications of adaptive PIC 
method, below we report two classic cases where uni- 
form PIC calculations show their limitations: coUi- 
sionless shocks and small scales objects (dust parti- 
cles) immersed in plasmas. 

5.1 Collisionless Shocks 

Simulations of collisionless shocks provide a sensitive 
test of the accuracy of particle rezoning methods [?] . 
In the slow shock calculations considered here, a mag- 
netized plasma is flowing toward a piston that reflects 
the particles. A switch off slow shock is considered, 
and the component of the magnetic field perpendic- 
ular to the normal of the piston is set to zero. 

We consider here the same conditions reported in 
Ref. [?] . The initial configuration is chosen according 
to the Rankinc-Hugoniot conditions. The initial ra- 
tios of the electron and ion pressures to the upstream 
magnetic field are (3e = A: = 0.01. The ratio of ion 
to electron mass is mi/rrie = 25; the ratio of the up- 
stream ion cyclotron and ion plasma frequencies is 
ujci/iLipi — 0.01, and the shock normal angle, with re- 
spect to the magnetic field, is V' = 75°. The size of 
the simulation region is L = 200 c/ujpi, and the shock 
is followed until ujdt = 50. Particles are injected at 
the right boundary to simulate a fiowing plasma. 

The simulations are performed using CE- 
LESTEID [?], a ID implicit PIC code, suitably 
modified by the author to include particle control. 



As a reference, we conduct a reference collision- 
less shock calculation with a uniform grid and with- 
out particle rezoning. The grid has 1000 cells giv- 
ing a uniform spacing with Aa; = 0.2 c/cUpi; 128 elec- 
trons and 128 ions per cell are used. Figure 4 shows 
the stack plot of as a function of the position at 
50 equally spaced time intervals between t — and 

LUcit = 50. 

The reference results are compared with a calcula- 
tion where particle rezoning is performed using the al- 
gorithms described above for splitting and for coales- 
cences. The computation uses an adaptive grid with 
finer spacing in the shock region (Ax ~ 0.5 c/wpi) 
and coarser outside. The region of fine spacing ex- 
pands in time to follow the motion of the shock. The 
grid spacing in the region of the shock is kept fixed; 
and, consequently, the grid spacing in the coarser re- 
gion grows to keep the number of grid points constant 
and equal to 300. Figure |4] shows the grid spacing at 
the end of the simulation, uucu = 50. Note that the 
area of the shock is well resolved, while the upstream 
region has large cells. We use the grid jiggling tech- 
nique of randomly displacing the grid spacing in the 
large cells to improve the energy conservation of the 
simulation [?]. This technique results in a random 
noise added to the grid spacing in the large cells. To 
avoid any noise in the shock region, the jiggling tech- 
nique is not used there. 

The particles are loaded with a uniform number per 
cell (the same as before), leading to higher accuracy 
where the grid is finer. Particle rezoning is required 
to keep the uniformity of the number of particles per 
cell as the grid is adapted. 

Figure [5] shows the profile of Bz at the end {uJcu = 
50) of the two calculation described above. Clearly, 
the evolution of the system is calculated correctly. In 
particular, the shock has traveled backward along the 
axis for a length of bOc/ojpi as in the reference case 
(Fig. 4) and as required by the Rankine-Hugoniot 
conditions. The results have been validated against 
previously published results [?] . 

5.2 Charging Of Dust Particles 

As a second test, we consider small objects (e.g. dust 
particles) immersed in a plasma. This condition is 
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Figure 4: Grid used in the adaptive shock calculation. 
The grid size in each cell is plotted versus the cell 
center. The region of smaller grid spacing moves to 
the left to follow the shock. 

common in industrial applications of plasma physics 
and in space and astrophysical occurrences of dusty 
plasmas. Dust particles immersed in plasmas tend 
to acquire a negative charge. The ions and electrons 
of the plasma reach the surface and stick to it. If 
no secondary emission or photoemission is present, 
the equilibrium charge on the dust particle must be 
negative to repel the more mobile electrons and at- 
tract the ions to achieve a balance of electron and 
ion currents. This problem is of interest in labo- 
ratory and in space plasmas [?, ?]. We consider 
here the case where a plasma with an ion to elec- 
tron temperature ratio T^jT^ — 20 and ion to elec- 
tron mass ratio mi/mf. — 1836 is drifting relative 
to a spherical dust particle of radius a/A/je = 0.4, 
where Xoe is the electron Debye length. The rel- 
ative velocity w is expressed by the Mach number 
M = wm^^"^ /{kTeY^^ = 10. The system is simu- 
lated using a cylindrical coordinate system with the 
vertical axis along the direction of the plasma flow 
and centered in the center of the spherical dust par- 
ticle. In this configuration, the azimuthal coordinate 
is invariant, and the problem is 2D axisymmetric. 
The interaction of the dust particle with the 
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Figure 5: Spatial profile of the z component of the 
magnetic field at the end of the simulation (nor- 
malized to its upstream value) at time ojdt = 50. 
Two different runs are shown. The first uses a uni- 
form 1000 cells grid (solid line) and the second an 
adaptive 300 cells grid (dashed line). 

plasma is described with the immersed boundary 
method. The application of the immersed boundary 
method in PIC codes is described in Ref. [?] for fiuid 
problems and in Ref. [?] for plasmas. In the present 
work, we will use the immersed boundary explicit 
PIC code DEMOCRITUS developed by the author 
for dusty plasma simulations [?] . A brief description 
of the method is given below, more details can be 
found in Ref. [?]. 

The dust particle is represented by motionless com- 
putational particles (object particles) with properties 
suitable to describe the macroscopic properties of the 
dust. Dust plasma interface conditions are treated 
with the immersed boundary method in two steps. 

First, we assign to the object particles a suscepti- 
bility Xp that can be interpolated to the vertices of 
the grid x^, to obtain a grid susceptibility: 

X„ = ^ S^pXp , (17) 
p 

where Syp are the linear assignment weights. The 
grid susceptibility is used to alter the Poisson's equa- 
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tion: 

T^cvil + Xv)Qvd(t)c' = Pc , (18) 

where the potential (j) and the charge density p are 
defined on the cell centers Xc and repeated indexes 
are summed. The operators Vcv and G^c are a differ- 
ence approximation of the divergence and gradient, 
respectively. As discussed in detail elsewhere [?, ?], 
Eq. (23) is solved everywhere, including in the inte- 
rior of the dust particle. The term (1 -I- Xv) gives 
an approximation to the correct interface conditions 
for the electric field. In the present case, Xv is the 
susceptibility of dielectric dust. 

Second, the object particles exert a friction on the 
plasma particles, via a slowing property fip that is 
interpolated to the grid, as in Eq. (22), to produce a 
grid quantity fiy used to introduce a damping term 
to the equation of motion of the plasma particles: 

= ^ EySyp — Vp ^ SypHv . (19) 
V V 

The second term in Eq. (24) can be as big as de- 
sired to stop the plasma particles on the surface of 
the dust. The damping term is zero everywhere out- 
side the region occupied by the dust. Equation (23) 
and Eq. (24) allow one to treat the field and particle 
boundary conditions on the surface of the dust. 

Figure [6] shows the configuration of the grid and 
of the dust particle for the problem considered here. 
Note that a nonuniform (but constant in time) grid 
is used to describe better the sheath around the dust 
particle. The distance of the dust particle from the 
boundaries is 10 Xoe- The plasma species are ini- 
tially loaded according to a drifting maxwellian dis- 
tribution with a downward vertical net flow velocity 
corresponding to a Mach number Af = 10. To reach 
an equilibrium, particles that flow out of the lower 
boundary are replaced by particles injected at the 
top boundary [?]. 

Figure [7] shows the history of the net charge accu- 
mulated on the dust particle. In this case, particle 
rezoning was used to ensure the accuracy of the cal- 
culation. The particles are loaded, initially, with a 
constant number of particles per cell, leading to a 
higher concentration around the dust particle where 
the cells are smaller. However, the plasma flow tends 
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Figure 6: Initial setup of a dust charging simulation. 
The dust particle is represented by material compu- 
tational particles with appropriate dielectric proper- 
ties for the immersed boundary method. An adaptive 
grid is used to resolve the small sub-Debye scale dust 
particle. 

to empty the region around the dust reducing the 
accuracy. Splitting the particles moving toward the 
dust and coalescing the particles moving away from 
it is desirable to keep the number of particles per cell 
and the accuracy constant. 

If the calculation is repeated without particle re- 
zoning, the accuracy worsens in time as the region 
around the dust becomes less populated. Two effects 
lead to decrease accuracy around the dust particle: 
the particles originally present are in part captured 
by the dust and in part just simply flow away ac- 
cording to their average downward velocity of Mach 
AI = 10. The new particles that replace them are 
flowing from regions of larger cells and are less nu- 
merous leading to a decrease of accuracy. 

As a result of the decrease in accuracy, the dust 
particle does not reach a steady state in the run with- 
out particle rezoning. 
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Figure 7: Evolution of the charge collected by the 
dust particle. Two runs are shown, both have uni- 
form grids but one has also particle control (solid 
line) and the other has no particle control (dashed). 



11 



